from __future__ import print_function

from Bio.Align.Applications import MuscleCommandline
from Bio import SeqIO
from Bio import AlignIO
from Bio.Seq import Seq
from Bio.SeqRecord import SeqRecord
from Bio.Alphabet import generic_dna
from Bio.Align import MultipleSeqAlignment
from Bio import pairwise2
from constructor.models import Node, Protocol, Primer
import logging
import os
import primers
import pydot
import re
import simplejson
import string
import subprocess
import sys
import time

logger = logging.getLogger('file_log')

"""This file should receive as input a target sequence, and output the steps
required to synthesize the sequence, as well as primers and target sequences
of each RPA reaction step."""

def createProtocol(paramsDict, seq_type, sequences, name, id):
	"""Main function."""
	namedTargets = sequences
	sequences = [target['sequence'] for target in sequences]
	target = [primers.parse_sequence(seq, seq_type) for seq in sequences]
	paramsDict['target'] = target
	fragments = paramsDict['fragmentList']

	logger.debug("\nNEW PROTOCOL REQUEST AT " + time.asctime(time.localtime()) + " with targets " + str(target) + "\n")
	
	if len(target) > 1:	
		targets = []
		for i in range(len(target)):
			targets.append(SeqRecord(target[i], id='sequence_' + str(i), description="sequence " + str(i)))
			
		child = subprocess.Popen(['muscle', '-clwstrict'], stdin=subprocess.PIPE, stdout=subprocess.PIPE, stderr=subprocess.PIPE, shell=False)
		SeqIO.write(targets, child.stdin, 'fasta')
		child.stdin.close()

		align = AlignIO.read(child.stdout, 'clustal')
		
		commonIndices, targets = get_common_indices(align, 50)

		logger.debug("\ngot common indices " + str(commonIndices) + ", targets= " + str(targets))

		treeList = []
		
		for sequence in targets:
			paramsDict['target'] = sequence
			
			if fragments[0] != 'false':
				matchThreshold = 30
				matchIndices = get_map_locations(sequence, paramsDict['fragmentList'], matchThreshold)
			
			else:
				matchThreshold = 30
				matchIndices = []
			
			matchIndices.extend(commonIndices[targets.index(sequence)])

			subSeqs = handle_fragments(sequence, paramsDict, matchThreshold, matchIndices)
			
			tree = construct_tree(subSeqs, paramsDict)
			tree = set_tree_depths(tree)
			
			treeList.append(tree)
			
		protocol = Protocol.objects.filter(id=id)[0]
		
		graph = pydot.Dot(graph_type='digraph')
		graph.set_label(name)
		
		oligos = []
		
		for tree in treeList:
			oligos.extend(get_oligos_from_tree(tree))
		
		for tree in treeList:
			graphResults = create_graph(tree, graph, protocol, namedTargets=namedTargets, combID=treeList.index(tree) + 1)
			graph = graphResults[0]
			protocol = graphResults[4]
			
		protocol.save()
		
		result = get_xdot(graph)
		
		protocol.graphString = result
		protocol.save()

		logger.debug("\n\n")
		
		return (result, protocol.id)
		
	else:
		paramsDict['target'] = target[0]
		target = target[0]

	if fragments[0] != 'false':
		matchThreshold = 30
		matchIndices = get_map_locations(target, paramsDict['fragmentList'], matchThreshold)
		subSeqs = handle_fragments(target, paramsDict, matchThreshold, matchIndices)
	
	else:
		subSeqs = generate_subsequences(paramsDict['target'], paramsDict['target'], 
										paramsDict['minOverlap'], paramsDict['maxOverlap'], 
										paramsDict['maxOligoLength'], paramsDict['maxDeviation'])
	
	tree = construct_tree(subSeqs, paramsDict)
	tree = set_tree_depths(tree)
	protocol = Protocol.objects.filter(id=id)[0]
	
	graph = pydot.Dot(graph_type='digraph')
	
	graphResults = create_graph(tree, graph, protocol, namedTargets=namedTargets)
	
	graph = graphResults[0]
	protocol = graphResults[4]
	protocol.save()
	graph.set_label(name)
		
	result = get_xdot(graph)
	
	protocol.graphString = result
	protocol.save()
	
	return (result, protocol.id)

def handle_fragments(target, paramsDict, matchThreshold, matchIndices):
	"""Returns the subsequences generated by the given set of natural fragments."""
	matchesWithOverlap = []
	for match in matchIndices:
		match = (match[0], match[1], 
				 get_best_fragment_overlap(paramsDict['target'], paramsDict['minOverlap'], 
										   paramsDict['maxOverlap'], match[0], match[1]), match[2])
		
		matchDict = {'sequence': target[match[0]:match[1]+1], 'startIndex': match[0], 
					 'endIndex': match[1], 'startOverlap': match[2][0], 'endOverlap': match[2][1], 
					 'nextEnd':'right', 'type': match[3]}
		
		if match[3] == 'combinatorial':
			bestFragDp, bestFragOverlap, bestFragScore = get_best_division(target[match[0]:match[1] + 1], paramsDict['minOverlap'],
							  paramsDict['maxOverlap'], paramsDict['maxOligoLength'], paramsDict['maxDeviation'], startIndex=match[0],)
			matchDict['bestFragDp'] = bestFragDp
			matchDict['bestFragOverlap'] = bestFragOverlap
		
		matchesWithOverlap.append(matchDict)

	# Sort the list of fragments by their start index, which allows us to split the list of fragments later on in generate_subsequences.
	sortedMatches = sorted(matchesWithOverlap, key=lambda k: k['startIndex'])

	logger.debug("starting with fragments: " + str(sortedMatches)) 
	
	return generate_subsequences(paramsDict['target'], paramsDict['target'], paramsDict['minOverlap'], paramsDict['maxOverlap'], 
								paramsDict['maxOligoLength'], paramsDict['maxDeviation'], fragments=sortedMatches)

def get_xdot(graph):
	process = subprocess.Popen(['dot', '-Txdot'], shell=False, stdin=subprocess.PIPE, stdout=subprocess.PIPE)
	process.stdin.write(graph.to_string())
	result, stderr = process.communicate()	
	return result
	
def get_oligos_from_tree(tree, remove=False):
	"""Returns a list of oligos in a given tree."""
	oligos = []
	
	if tree['leftRxn'] and tree['rightRxn']:
		oligos.extend(get_oligos_from_tree(tree['leftRxn']))
		oligos.extend(get_oligos_from_tree(tree['rightRxn']))
		return oligos
		
	else:
		return [tree['target']]

def set_tree_depths(tree, depth=0):
	"""Adds a key-value for tree depth to the tree dict."""
	tree['depth'] = depth
	
	if tree['leftRxn']:
		tree['leftRxn'] = set_tree_depths(tree['leftRxn'], depth + 1)
	
	if tree['rightRxn']:
		tree['rightRxn'] = set_tree_depths(tree['rightRxn'], depth + 1)
		
	return tree
	
def get_number_nodes(tree, count=0):
	"""Gets the total number of nodes in a tree."""
	count += 1
	
	if tree['leftRxn'] and tree['rightRxn']:
		count += get_number_nodes(tree['leftRxn'])
		count += get_number_nodes(tree['rightRxn'])
		
	return count
	
def construct_tree(subSeqs, paramsDict, parentNode=False, side=False, startIndex=0):
	"""Builds a tree (in the form of a dict) out of a given set of subsequences
	and optimizes primer size for each reaction."""
	node = {'target':'', 'type':'', 'side':'', 'forwardPrimer':'', 'reversePrimer':'', 'startIndex':''}
	subSeq = subSeqs[0]
	node['target'] = str(subSeq['sequence'])
	node['type'] = subSeq['type']
	node['side'] = side
	
	if not(side) or side == 'left':
		if not(parentNode):
			node['target'] = str(subSeq['sequence'])
		
		else:
			node['target'] = str(subSeq['sequence'])
		
	else:
		node['target'] = str(subSeq['sequence'].reverse_complement())
   
	# If parent node does not exist, calculate both forward and reverse primers.
	if not(parentNode):
		primerResults = primers.main(paramsDict, 'string', node['target'])
		
		forwardPrimerSize = primerResults[0]
		reversePrimerSize = primerResults[2]
		
		node['forwardPrimer'] = Primer(sequence=node['target'][0:forwardPrimerSize], size=forwardPrimerSize)
		node['forwardPrimer'].save()

		node['reversePrimer'] = Primer(sequence=str(Seq(node['target'][len(node['target']) - reversePrimerSize:]).reverse_complement()),
									   size=reversePrimerSize)
		node['reversePrimer'].save()
		
		node['startIndex'] = startIndex
	
	# If parent node exists, inherit the appropriate primer from the parent node and calculate the other.
	else:
		if side == 'right':
			offset = parentNode['startIndex'] - 1 if parentNode['startIndex'] != 0 else 0
			node['startIndex'] = parentNode['dp'] - parentNode['overlap'] / 2 + 1 + offset
			if len(subSeqs) > 1:
				node['reversePrimer'] = parentNode['reversePrimer']
				
				primerResults = primers.main(paramsDict, 'string', node['target'], forwardOnly=True)
				
				forwardPrimerSize = primerResults[0]
				node['forwardPrimer'] = Primer(sequence=str(Seq(node['target'][len(node['target']) - forwardPrimerSize:]).reverse_complement()),
											   size=forwardPrimerSize)
				node['forwardPrimer'].save()
			
		elif side == 'left':
			node['startIndex'] = parentNode['startIndex']
			if len(subSeqs) > 1:
				node['forwardPrimer'] = parentNode['forwardPrimer']
				
				primerResults = primers.main(paramsDict, 'string', node['target'], reverseOnly=True)
				
				reversePrimerSize = primerResults[0]
				node['reversePrimer'] = Primer(sequence=str(Seq(node['target'][len(node['target']) - reversePrimerSize:]).reverse_complement()),
											   size=reversePrimerSize)
				node['reversePrimer'].save()
	
	# If we're not at an oligo, store dp and overlap and continue generating the tree.
	if len(subSeqs) > 1:
		sequences = [dict['sequence'] for dict in subSeqs]
		sequenceEnds = [dict['endIndex'] for dict in subSeqs]
		rightSeqIndex = sequenceEnds[1:].index(sequenceEnds[0]) + 1
		node['dp'] = subSeq['bestDp']
		node['overlap'] = subSeq['bestOverlap']
		
		node['leftRxn'] = construct_tree(subSeqs[1:rightSeqIndex], paramsDict, node, 'left')
		node['rightRxn'] = construct_tree(subSeqs[rightSeqIndex:], paramsDict, node, 'right')
		
	else:
		node['leftRxn'] = False
		node['rightRxn'] = False
		node['forwardPrimer'] = False
		node['reversePrimer'] = False
	 
	return node

def hide_nodes(dotData, name):
	dotData = dotData.replace("\t", "")

	print(dotData, file=sys.stderr)

	graph = pydot.graph_from_dot_data(dotData)
	print('edges: ', len(graph.get_edges()), file=sys.stderr)
	node = graph.get_node(name)[0]
	childNames = get_child_names(graph, name)
	childNames.append(name)

	print('nodes: ', len(graph.get_nodes()), file=sys.stderr)

	for edge in graph.get_edges():
		print('edge source: ', edge.get_source(), file=sys.stderr)

	for childName in childNames:
		child = graph.get_node(childName)[0]
		if child.get('fontcolor') == 'white' and node.get('color') == 'white':
			child.set('fontcolor', '')
			child.set('color', child.get('comment'))
			child.set('href','javascript:void(0)')
		
		else:
			child.set('comment', child.get('color'))
			child.set('fontcolor', 'white')
			child.set('color', 'white')
			child.set('href', '')
			
	for edgeTuple in get_child_edges(graph, childNames):
		edge = graph.get_edge(edgeTuple[0], edgeTuple[1])[0]

		print >> sys.stderr, 'edge:'
		print >> sys.stderr, edge

		if edge.get('fontcolor') == 'white' and edge.get('color') == 'white':
			edge.set('fontcolor', 'grey')
			edge.set('color', edge.get('comment'))
			edge.set('href','javascript:void(0)')
		
		else:
			edge.set('comment', edge.get('color'))
			edge.set('fontcolor', 'white')
			edge.set('color','white')
			edge.set('href','')
	
	return get_xdot(graph)

def get_child_names(graph, name):
	"""Gets the names of all descendant nodes of a given node."""
	children = []
	for edge in graph.get_edges():
		if edge.get_source() == name:
			destination = graph.get_node(edge.get_destination())[0]
			children.append(destination.get_name())
			children.extend(get_child_names(graph, destination.get_name()))
		
	return children
	
def get_child_edges(graph, childNames):
	"""Gets all child edges of a node, given a list of its child nodes."""
	childEdges = []
	for edge in graph.get_edges():
		if edge.get_source() in childNames and edge.get_destination() in childNames:
			childEdges.append((edge.get_source(), edge.get_destination()))
			
	return childEdges

def create_graph(tree, graph, protocol, parentNode=False, oligos=0, count=0, combID=False, namedTargets=False):
	"""Creates a graph object and saves it in the database."""
	nodeName = str(tree['target'])
	
	if tree['type'] != 'common' and (tree['leftRxn'] and tree['rightRxn']):
		nodeName += str(combID)
	
	if len(graph.get_node(nodeName)) > 0 and not(tree['leftRxn'] and tree['rightRxn']) and tree['type'] != 'common':
		tree['type'] = 'subCommon'
		
	node = pydot.Node(nodeName) if len(graph.get_node(nodeName)) == 0 else graph.get_node(nodeName)[0]

	node.set('style', 'bold')
	count += 1
	node.set('shape', 'box')
	node.set('href', 'javascript:void(0)')
	
	if combID:
		prevIDs = node.get('id') + ', V' + str(combID) if node.get('id') else 'V' + str(combID)
		node.set('id', prevIDs)
		if tree['type'] != 'common' and tree['type'] != 'subCommon':
			node.set('colorscheme', 'set19')
			node.set('color', str(combID))
			node.set('comment', node.get('color'))
		else:
			node.set('color', 'black')
			node.set('comment', 'black')
	
	if len(graph.get_node(nodeName)) == 0:
		graph.add_node(node)
	
	nodeDB = Node()
	nodeDB.protocol = protocol
	nodeDB.target = tree['target']
	nodeDB.treeDepth = tree['depth']
	
	if tree['forwardPrimer'] and tree['reversePrimer']:
		nodeDB.forwardPrimer = tree['forwardPrimer']
		nodeDB.reversePrimer = tree['reversePrimer']
		
	nodeDB.startIndex = tree['startIndex']
	nodeDB.isNaturalFragment = True if tree['type'] == 'fragment' else False
	nodeDB.isRoot = False if parentNode else True
	nodeDB.isCommon = True if tree['type'] == 'common' or tree['type'] == 'subCommon' else False

	if not(tree['leftRxn']) and not(tree['rightRxn']):
		oligos += 1
		
		label = 'Oligo ' + str(oligos)
		
		if nodeDB.isNaturalFragment:
			label += ' (NF)'
			
		if tree['type'] == 'common':
			label = ''.join(['Common ', label])
		
		if combID and not(nodeDB.isCommon):
			label += ' V' + str(combID)
			
		if tree['type'] == 'subCommon':
			label += ' ('
			for sameNode in graph.get_node(nodeName):
				label += sameNode.get_id()
			
			label += ')'

		node.set_label(label)
			 
		nodeDB.leftReaction = None
		nodeDB.rightReaction = None
		nodeDB.isLeaf = True
		
	else:
		nodeDB.divPoint = tree['dp']
		(graph, oligos, leftNode, nodeDB.leftReaction, protocol, count) = (
							create_graph(tree['leftRxn'], graph, protocol, node, oligos, count, 
										 combID=combID))
		
		(graph, oligos, rightNode, nodeDB.rightReaction, protocol, count) = (
							create_graph(tree['rightRxn'], graph, protocol, node, 
										 oligos, count, combID=combID))
		
		label = (get_label(leftNode.get('label'), rightNode.get('label')) if parentNode else 'Target Sequence')
		
		if nodeDB.isNaturalFragment:
			label += ' (NF)'
		
		if combID and not(nodeDB.isCommon):
			label += ' V' + str(combID)
			
		if tree['type'] == 'common':
			label = ''.join(['Common ', label])
			
		elif tree['type'] == 'subCommon':
			label += ' ('
			for sameNode in graph.get_node(nodeName):
				label += sameNode.get('id')
			
			label += ')'
			
		node.set_label(label)
		
	nodeDB.name = node.get_label().replace(' ', '_')
	nodeDB.protocol = protocol
	nodeDB.save()

	if parentNode:
		if tree['side'] == 'right' and tree['leftRxn']:
			nodeDB.forwardPrimer.name = nodeDB.name + '_Forward_Primer'
			nodeDB.forwardPrimer.save()
		
		elif tree['side'] == 'left' and tree['leftRxn']:
			nodeDB.reversePrimer.name = nodeDB.name + '_Reverse_Primer'
			nodeDB.reversePrimer.save()
		
		if len(graph.get_edge(parentNode.get_name(), nodeName)) == 0:
			edge = pydot.Edge(parentNode, node)
			edge.set('href', 'javascript:void(0)')
			
			onclick = """\"edgeClick('%s')\""""
			edge.set('onclick', onclick % ''.join(node.get_name().replace('"', '').split()))
			edge.set_label('+')
			edge.set_fontcolor('grey')
			graph.add_edge(edge)
	
	else:
		nodeDB.forwardPrimer.name = nodeDB.name + '_Forward_Primer'
		nodeDB.reversePrimer.name = nodeDB.name + '_Reverse_Primer'
		nodeDB.forwardPrimer.save()
		nodeDB.reversePrimer.save()

		targetIndex = [target['sequence'] for target in namedTargets].index(tree['target'])

		if namedTargets[targetIndex]['name'] != 'NOTAVARIABLE':
			node.set_label(namedTargets[targetIndex]['name'])

		protocol.rootNode = nodeDB
		
	protocol.save()
	node.set('onclick', '"nodeClick(' + str(nodeDB.id) + ')"')
	return (graph, oligos, node, nodeDB, protocol, count)

def get_label(left, right):
	"""Used by create_graph to get the properly formatted name of a node given its left and right children."""
	leftOligo = re.findall(r"\d+", left)[0]
	rightOligo = re.findall(r"\d+", right)
	
	rightOligo = rightOligo[0 if right.find('Oligos') == -1 else 1]
	
	return 'Oligos ' + leftOligo + '-' + rightOligo

def get_best_fragment_overlap(target, minOverlap, maxOverlap, startIndex, endIndex, verbose=False):
	"""Finds the best overlap between a fragment and its adjacent oligos given where it starts and ends on target."""
	count = 1
	
	for overlapLength in range(minOverlap, maxOverlap + 1):
		if startIndex != 0: # No point trying to calculate overlap in the beginning if our fragment is at the start of target.
			leftOverlap = (Seq((str(target)[startIndex:startIndex + overlapLength]), generic_dna))
			
			# Calculate the total score for the overlap between left end of fragment and rest of target.
			scoreLeftF = primers.get_score_from_sequences(leftOverlap, 
						target[0:startIndex - 1], oneAlignment=True, debug=verbose)
			
			scoreRightF = primers.get_score_from_sequences(leftOverlap, 
						target[startIndex + overlapLength:], oneAlignment=True, debug=verbose)
			
			scoreLeftR = primers.get_score_from_sequences(leftOverlap, 
						target.complement()[0:startIndex - 1], debug=verbose, oneAlignment=True)
			
			scoreRightR = primers.get_score_from_sequences(leftOverlap, 
						target.complement()[startIndex + overlapLength:], debug=verbose, oneAlignment=True)
			
			selfScore = primers.get_score_from_sequences(leftOverlap, leftOverlap) # Score of desired reaction.
			
			leftValid = True
			
			if max(scoreLeftF, scoreRightF, scoreLeftR, scoreRightR) >= selfScore:
				leftValid = False
			
			leftOverlapScore = scoreLeftF + scoreRightF + scoreLeftR + scoreRightR
			
			if count == 1 or (leftOverlapScore < bestLeftScore and leftValid): # Save the current score if it is better than the previous best.
				bestLeftScore = leftOverlapScore
				bestLeftOverlap = overlapLength
				bestLeftValid = leftValid
				
		else:
			bestLeftOverlap = False

		if endIndex != len(target) - 1: # No point calculating right overlap if fragment goes to the end of target.
			rightOverlap = (Seq((str(target)[endIndex - overlapLength:endIndex]), generic_dna))
			
			# Calculate the total score for the overlap between left end of fragment and rest of target.
			scoreLeftF = primers.get_score_from_sequences(rightOverlap, 
						target[0:endIndex - overlapLength], oneAlignment=True, debug=verbose)
			
			scoreRightF = primers.get_score_from_sequences(rightOverlap, 
						target[endIndex + 1:], oneAlignment=True, debug=verbose)
			
			scoreLeftR = primers.get_score_from_sequences(rightOverlap, 
						target.complement()[0:endIndex - overlapLength], debug=verbose, oneAlignment=True)
			
			scoreRightR = primers.get_score_from_sequences(rightOverlap, 
						target.complement()[endIndex + 1:], debug=verbose, oneAlignment=True)
			
			selfScore = primers.get_score_from_sequences(rightOverlap, rightOverlap) # Score of desired reaction.
			
			rightValid = True
			
			if max(scoreLeftF, scoreRightF, scoreLeftR, scoreRightR) >= selfScore:
				rightValid = False
			
			rightOverlapScore = scoreLeftF + scoreRightF + scoreLeftR + scoreRightR
			
			if count == 1 or (rightOverlapScore < bestRightScore and rightValid): # Save the current score if it is better than the previous best.
				bestRightScore = rightOverlapScore
				bestRightOverlap = overlapLength
				bestRightValid = rightValid
				
		else:
			bestRightOverlap = False

		count += 1
		
	return bestLeftOverlap, bestRightOverlap

def get_map_locations(target, fragmentList, matchThreshold):
	"""Returns a list of starting and ending indices where each fragment matches
	target at greater than matchThreshold bases in a row."""
	for fragment in fragmentList:
		alignment = pairwise2.align.localms(target, fragment, 2, -1, -.5, -.1, one_alignment_only=True)
		
		for align1, align2, score, begin, end in alignment:
			i = 0
			while align1[i] == '-':
				i += 1
			indices = primers.get_indices_from_alignment(align1[i:], align2[i:], debug=True)[0]
		
		matchIndices = get_start_and_end_indices(indices, matchThreshold)
		matchIndices = [(match[0], match[1], 'natural') for match in matchIndices]
	
	return matchIndices

def get_start_and_end_indices(indices, threshold):
	"""Given a list of indices where two sequences have
	the same nucleotides, returns the start and end indices
	for every matching region of more than threshold bases."""
	start = 0
	matchesInRow = 0
	matchIndices = [] # matchIndices will be a list of tuples, where tuple[0] is the fragment, tuple[1] is the start of a match, tuple[2] is the end
	for i in range(len(indices) - 1):
		if (int)(indices[i + 1]) == (int)(indices[i]) + 1:
			if matchesInRow == 0:
				start = indices[i]
			last = indices[i + 1]
			matchesInRow += 1
			
		else:
			if matchesInRow > threshold:
				matchIndices.append((start, last))
			matchesInRow = 0
			
	if matchesInRow > threshold:
		matchIndices.append((start, last))
	
	return matchIndices

def get_common_indices(align, matchThreshold):
	"""Given a multiple sequence alignment object, returns a list
	of areas where all targets have the same sequence in the form
	of a list of start indices and end indices."""
	sequences = [record.seq for record in align]
	
	# Clustal can reorder target sequences, so we need to return a list of target sequences in the same order as they are in the alignment.
	allchars = ''.join([chr(x) for x in range(256)])
	targetsReordered = [Seq(str(seq).translate(allchars, '-'), generic_dna) for seq in sequences]
	matchLocations = []
	for i in range(len(sequences[0])): # iterate through each sequence, finding the nucleotide indices they all have in common
		bases = [seq[i] for seq in sequences]
		if bases[1:] == bases[:-1]:
			matchLocations.append(i)
		
	matchLocations = get_start_and_end_indices(matchLocations, matchThreshold)
	adjustedMatches = []
	
	for i in range(len(sequences)): # adjust for the presence of dashes in the alignment string, causing each sequence to have slightly different indices
		seq = sequences[i]
		adjustedMatches.append([])
		for j in range(len(matchLocations)):
			match = matchLocations[j]
			newMatch = (match[0] - str(seq[0:match[0]]).count('-'), 
						match[1] - str(seq[0:match[1]]).count('-'),
						'combinatorial')
			adjustedMatches[i].append(newMatch)
	return (adjustedMatches, targetsReordered)

def get_best_division(target, minOverlap, maxOverlap, maxOligoLength, maxDeviation, verbose=False, startIndex=0):
	"""Finds the best division point, starting in the center of target and working outwards."""
	subSeqs = []
	dp = len(target) / 2
	count = 1
	minScore = 0
	
	while len(target) / 2 - maxDeviation < dp < len(target) / 2 + maxDeviation:
		for overlapLength in range(minOverlap, maxOverlap + 1):
			
			leftSeq = Seq(str(target)[0:dp + (overlapLength + 1) / 2], generic_dna)
			rightSeq = Seq(str(target)[dp - overlapLength / 2:], generic_dna)
			
			overlap = (Seq((str(target)[dp - overlapLength / 2:
										dp + (overlapLength + 1) / 2]), generic_dna))
			
			# Calculate the total score of potential non-specific interactions between sequences.
			scoreLeftF = primers.get_score_from_sequences(overlap, 
						target[0:dp - overlapLength / 2], oneAlignment=True, debug=verbose)
			
			scoreRightF = primers.get_score_from_sequences(overlap, 
						target[dp + overlapLength / 2:], oneAlignment=True, debug=verbose)
			
			scoreLeftR = primers.get_score_from_sequences(overlap, 
						target.complement()[0:dp - overlapLength / 2], debug=verbose, oneAlignment=True)
			
			scoreRightR = primers.get_score_from_sequences(overlap, 
						target.complement()[dp + overlapLength / 2:], debug=verbose, oneAlignment=True)
			
			selfScore = primers.get_score_from_sequences(overlap, overlap) # Score of desired reaction.
			
			valid = True
			
			if max(scoreLeftF, scoreRightF, scoreLeftR, scoreRightR) >= selfScore:
				valid = False
			
			score = scoreLeftF + scoreRightF + scoreLeftR + scoreRightR
			
			if verbose:
				print("\n" + leftSeq + "\n" + (" ")*(dp - overlapLength / 2) + overlap + "\n" + (" ")*(dp - overlapLength / 2) + rightSeq + "\n")
				print("dp: " + str(dp) + " overlapLength: " + str(overlapLength) + " count: " + str(count) + "\noverlap: " + overlap + 
					  "\nhas score " + str(score) + "\n\n")
				
			if count == 1 or (score < bestScore and valid): # Save the current score if it is better than the previous best.
				bestScore = score
				bestDp = dp
				bestOverlap = overlapLength
				bestValid = valid
			
		dp += count * int(pow(-1, count % 2)) # DP will go back and forth from the starting position.
		count += 1
		
		return bestDp, bestOverlap, bestScore

def generate_subsequences(fullTarget, target, minOverlap, maxOverlap, 
						  maxOligoLength, maxDeviation, verbose=False, 
						  startIndex=0, fragments=False, isCommon=False):
	"""Generates a list of sub-sequences to be synthesized."""
	
	if len(target) == 0:
		return []

	subSeqs = []
	dp = len(target) / 2
	count = 1
	minScore = 0

	logger.debug("startindex: " + str(startIndex) + ", endindex: " + str(startIndex - 1 + len(str(target))) + ", fragments: " + str(fragments))

	if len(target) < maxOligoLength:
		if not(fragments) and not(isCommon):
			type = 'oligo'
		
		elif fragments:
			type = 'fragment'
		
		else:
			type = 'common'
			
		return [{'sequence': target, 'startIndex': startIndex, 'bestDp': '', 'bestOverlap': '', 
				 'endIndex': startIndex + len(target) - 1, 'type': type}]
	
	else:
		bestDp, bestOverlap, bestScore = get_best_division(target, minOverlap, maxOverlap, 
								maxOligoLength, maxDeviation, verbose=False, startIndex=startIndex)
		if fragments:
			splitIndex = False
			lastEnd = False
			lastIndexBeforeSplit = False
			for fragment in fragments:
				if fragment['startIndex'] < bestDp + startIndex < fragment['endIndex'] and fragment['endIndex'] > lastEnd:
					splitIndex = fragments.index(fragment)
					lastEnd = fragments[splitIndex]['endIndex']
				
				elif fragment['endIndex'] < bestDp + startIndex:
					lastIndexBeforeSplit = fragments.index(fragment)
			
			logger.debug("splitindex: " + str(splitIndex) + ", bestDp: " + str(bestDp))

			if splitIndex != False or splitIndex == 0:
				fragment = fragments[splitIndex]
				if fragment['nextEnd'] == 'right':
					fragment['nextEnd'] = 'left'
					subSeqs.append({'sequence': target, 'startIndex': startIndex, 
									'bestDp': fragment['endIndex'] - fragment['endOverlap'],
									'bestOverlap': fragment['endOverlap'], 'endIndex': startIndex + len(target) - 1, 'type': 'oligo'})
					
					containedFragments = []

					if fragment['type'] == 'combinatorial':
						for contFrag in fragments:
							if fragment['startIndex'] < contFrag['startIndex'] < fragment['endIndex'] \
								and contFrag['endIndex'] < fragment['endIndex']:
								containedFragments.append(contFrag)

					subSeqs.extend(generate_subsequences(fullTarget, Seq(str(target)[0:fragment['endIndex'] - startIndex + 1], generic_dna),
														minOverlap, maxOverlap, maxOligoLength, maxDeviation,
														verbose, startIndex, fragments[0:splitIndex + 1]+ containedFragments))
					
					# Generate subsequences from the end of the current fragment (plus overlap) to the end of the target.
					subSeqs.extend(generate_subsequences(fullTarget, Seq(str(target)[fragment['endIndex'] - startIndex + 1 - fragment['endOverlap']:], generic_dna), 
													 minOverlap, maxOverlap, maxOligoLength, maxDeviation, verbose, 
													 fragment['endIndex'] - fragment['endOverlap'] + 1, 
													 fragments[splitIndex + 1:] if len(fragments) > 0 else False))
				
				# If we are at the left end of a fragment, time to divide up everything to the right of that end.
				else:				   
					if fragment['endIndex'] != len(fullTarget) - 1:
						if str(target) != str(fragment['sequence']):
							type = 'oligo'
						elif fragment['type'] == 'combinatorial':
							type = 'common'
						else:
							type = 'fragment'
						subSeqs.append({'sequence': target, 'startIndex': startIndex, 
										'bestDp': fragment['startIndex'] - startIndex + 1 + fragment['startOverlap'], 
										'bestOverlap': fragment['startOverlap'], 
										'endIndex': fragment['endIndex'], 'type': type})
					
					subSeqs.extend(generate_subsequences(fullTarget, Seq(str(target)[0:fragment['startIndex'] - startIndex + fragment['startOverlap']], generic_dna), 
													 minOverlap, maxOverlap, maxOligoLength, maxDeviation, 
													 verbose, startIndex, fragments[0:splitIndex] if len(fragments) > 0 else False))
					
					if fragment['type'] == 'natural': # If the current fragment is a natural fragment, don't split it.
						if str(target) != str(fragment['sequence']):
							subSeqs.append({'sequence': fragment['sequence'], 'startIndex': fragment['startIndex'], 
											'bestDp': fragment['startIndex'] + fragment['startOverlap'] / 2, 
											'bestOverlap': fragment['startOverlap'], 'endIndex': fragment['endIndex'], 'type': 'fragment'})
					
					else: # If the fragment is a common sequence in a combinatorial protocol, divide it up like a normal sequence.
						bestFragDp = fragment['bestFragDp']
						bestFragOverlap = fragment['bestFragOverlap']
						if str(target) != str(fragment['sequence']): # POTENTIAL BUG when fragment is equal to target
							# Check to see if there are any fragments contained in the current one.
							containedFragments = []
							for frag in fragments:
								if frag['startIndex'] > fragment['startIndex'] and frag['endIndex'] < fragment['endIndex']:
									containedFragments.append(frag)

							if len(containedFragments) == 0:
								subSeqs.append({'sequence': fragment['sequence'], 'startIndex': fragment['startIndex'],
											  'bestDp': fragment['bestFragDp'], 'bestOverlap': fragment['bestFragOverlap'], 
											  'endIndex': fragment['endIndex'], 'type': 'oligo'})
							
								subSeqs.extend(generate_subsequences(fullTarget, Seq(str(fragment['sequence'])[0:bestFragDp + (bestFragOverlap + 1) / 2], generic_dna), minOverlap,
													 maxOverlap, maxOligoLength, maxDeviation, verbose, fragment['startIndex'], isCommon=True))
								
								subSeqs.extend(generate_subsequences(fullTarget, Seq(str(fragment['sequence'])[bestFragDp - bestFragOverlap / 2:], generic_dna), minOverlap,
													 maxOverlap, maxOligoLength, maxDeviation, verbose, fragment['startIndex'] + bestFragDp - bestFragOverlap / 2, isCommon=True))

							else:
								subSeqs.extend(generate_subsequences(fullTarget, fragment['sequence'], minOverlap, 
													maxOverlap, maxOligoLength, maxDeviation, verbose, 
													fragment['startIndex'], containedFragments, isCommon=True))

						else:
							if len(str(fragment['sequence'])) > maxOligoLength:
								subSeqs.extend(generate_subsequences(fullTarget, Seq(str(fragment['sequence'])[0:bestDp + (bestOverlap + 1) / 2], generic_dna), minOverlap,
																		 maxOverlap, maxOligoLength, maxDeviation, verbose, fragment['startIndex'], isCommon=True))
									
								subSeqs.extend(generate_subsequences(fullTarget, Seq(str(fragment['sequence'])[bestDp - bestOverlap / 2:], generic_dna), minOverlap,
																	 maxOverlap, maxOligoLength, maxDeviation, verbose, fragment['startIndex'] + bestDp - bestOverlap / 2, isCommon=True))
					
			elif fragments[0]['endIndex'] - startIndex < bestDp:			
				subSeqs.append({'sequence': target, 'startIndex': startIndex, 'bestDp': bestDp, 'bestOverlap': bestOverlap, 
								'endIndex': startIndex + len(target) - 1, 'type': 'oligo'})
				
				subSeqs.extend(generate_subsequences(fullTarget, Seq(str(target)[0:bestDp + (bestOverlap + 1) / 2], generic_dna), minOverlap,
												 maxOverlap, maxOligoLength, maxDeviation, verbose, startIndex, [fragments[0]]))
				
				subSeqs.extend(generate_subsequences(fullTarget, Seq(str(target)[bestDp - bestOverlap / 2:], generic_dna), minOverlap,
												 maxOverlap, maxOligoLength, maxDeviation, verbose, 
												 startIndex + bestDp - bestOverlap / 2, fragments[1:] if len(fragments) > 0 else False))
				
			else:
				subSeqs.append({'sequence': target, 'startIndex': startIndex, 'bestDp': bestDp, 'bestOverlap': bestOverlap, 
								'endIndex': startIndex + len(target) - 1, 'type': 'oligo'})
				
				subSeqs.extend(generate_subsequences(fullTarget, Seq(str(target)[0:bestDp + (bestOverlap + 1) / 2], generic_dna), minOverlap,
												 maxOverlap, maxOligoLength, maxDeviation, verbose, startIndex))
				
				subSeqs.extend(generate_subsequences(fullTarget, Seq(str(target)[bestDp - bestOverlap / 2:], generic_dna), 
													 minOverlap, maxOverlap, maxOligoLength, maxDeviation, verbose, 
													 startIndex + bestDp - bestOverlap / 2, fragments))
				
		else:			
			subSeqs.append({'sequence': target, 'startIndex': startIndex, 'bestDp': bestDp, 'bestOverlap':bestOverlap, 
							'endIndex': startIndex + len(target) - 1, 'type': 'oligo' if not(isCommon) else 'common'})
			
			subSeqs.extend(generate_subsequences(fullTarget, Seq(str(target)[0:bestDp + (bestOverlap + 1) / 2], generic_dna), minOverlap,
											 maxOverlap, maxOligoLength, maxDeviation, verbose, startIndex, isCommon=isCommon))
			
			subSeqs.extend(generate_subsequences(fullTarget, Seq(str(target)[bestDp - bestOverlap / 2:], generic_dna), minOverlap,
											 maxOverlap, maxOligoLength, maxDeviation, verbose, startIndex + bestDp - bestOverlap / 2, isCommon=isCommon))
		
		return subSeqs
